---
title: "CARRS RDD analysis using a standardized cutoff"
knit: (function(input, ...) {rmarkdown::render(input, output_file=".pdf")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"



---
```{r , setup, include=FALSE, eval = TRUE}

knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""

)


 knitr::opts_knit$set(
   root.dir = ""    
 )
```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
  library(ggpubr)
```

```{r loaddata}
data<-data.frame(read.csv("carrs_bp.csv"))
```

```{r fuc}
##Defining the common function##
  carrs_rdd = function(a,b,d,outcome,od,yl,boxh,boxv,clr,l,u) {
    
    d$Y<-a
    d$X<-b
    gname=od
    yl=yl
    l=l
    u=u
    runningvar=b
    out<- rdrobust(d$Y, d$X, c = 0, 
                    cluster =d$pid, 
                    covs = cbind(d$w02, d$w24), 
                    kernel = "triangular", 
                    p = 1,rho=1)

    ret <- data.frame(
    sex = c(formatC(round(out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(out$ci[3,1],2),2,format="f"), " - ", formatC(round(out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2])))
    
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
    output<<-ret %>% add_rownames 


    d$bw <- out$bws[1,1]
    bw = out$bws[1,1]

   plotsub<-d[runningvar>=-bw & runningvar<=bw,] ##the plot is only for the values within the bandwidth
  

   c=0
  
  
  plot1 =rdplot(plotsub$Y, plotsub$X,
                  kernel = "triangular", p=1, 
                  covs = cbind(plotsub$w02, plotsub$w24), 
                         ci=95, hide=TRUE)
  rdplot_mean_bin = plot1$vars_bins[,"rdplot_mean_bin"]
  rdplot_mean_y   = plot1$vars_bins[,"rdplot_mean_y"]
  y_hat           = plot1$vars_poly[,"rdplot_y"]
  x_plot          = plot1$vars_poly[,"rdplot_x"]
  rdplot_cil_bin =  plot1$vars_bins[,"rdplot_ci_l"]
  rdplot_cir_bin =  plot1$vars_bins[,"rdplot_ci_r"]
  rdplot_mean_bin=  plot1$vars_bins[,"rdplot_mean_bin"]
  
  rdplot_mean_y = rdplot_mean_y
  y_hat = y_hat

  y_hat_r=y_hat[x_plot>=c & x_plot!=0]
  y_hat_l=y_hat[x_plot<c & x_plot!=0]
  x_plot_r=x_plot[x_plot>=c & x_plot!=0]
  x_plot_l=x_plot[x_plot<c & x_plot<0]
  

  x.label="\nStandardized blood pressure\n"
  y.label=outcome
  x.lim=c(min(plotsub$X, na.rm=T),max(plotsub$X, na.rm=T))
  y.lim=c(l,u)
  
  
  
  color_plot<- ggplot() + 
    geom_point(aes(x = rdplot_mean_bin, y = rdplot_mean_y, 
                    ), size=2,  colour ="#999999", na.rm = TRUE) +
    geom_line(aes(x = x_plot_l, y = y_hat_l,colour ='x_plot_l'),
                col= clr, size=1,na.rm = TRUE) +
    geom_line(aes(x = x_plot_r, y = y_hat_r,colour ='x_plot_r'), 
                col= clr, size=1, na.rm = TRUE) +
    
    scale_x_continuous(breaks=seq(-3,3,0.4)) +
    labs(x = x.label, y = y.label) +
    labs(caption = paste("MSE optimal bandwidth: +/- ", round(bw,2),"\n")) +
    theme(legend.position = "None") +
    geom_vline(xintercept = c, size = 0.5) +
    
    coord_cartesian(xlim = x.lim, ylim = y.lim) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          axis.text=element_text(family="serif", size=28), 
          axis.title=element_text(family="serif", size=28),
          axis.title.x = element_text(vjust = 5),
          axis.title.y = element_text(vjust = -5),
          plot.caption = element_text(family="serif", size=24, vjust = 10))


    out<<-out
    bw<<-bw 
    color_plot<<- color_plot
    plot1<<-plot1
    y_hat_l<<-y_hat_l
    y_hat_r<<-y_hat_r
    x_plot_l<<-x_plot_l
    x_plot_r<<-x_plot_r
    plotsub<<-plotsub
  }

```

**Graphs

\newpage

**Outcome : Systolic**

**Females and males**  



```{r dfull}
data<-as.data.table(data)

carrs_rdd(data$sysdiff,data$runningvar,data,"Change (mmHg)\n\n","overall_bp_abs.pdf",0.1,-0.9,10.5,"#25AC82FF",-25,25)
sysdiff_total <-output%>%mutate(Total=sex)%>%select(-sex)
bwTsysdiff  <- bw 
plotsysdiffT<-color_plot
rm(bw)

```
  


```{r ,results= TRUE}
summary(out)
```
\newpage

**Outcome : Systolic**

**Female**  


```{r dfemale}
data1<-data[female==1,]

carrs_rdd(data1$sysdiff,data1$runningvar_female,data1,"Change (mmHg)\n\n","female_bp_abs.pdf",0.1,-0.85,10.5,"#46337EFF",-25,25)
sysdiff_female <-output%>%mutate(Female=sex)%>%select(-sex)
bwFsysdiff  <- bw
plotsysdiffF<-color_plot
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```
\newpage

**Outcome : Systolic (Males)**

```{r dmale}
data2<-data[female==0,]

carrs_rdd(data2$sysdiff,data2$runningvar_male,data2,"Change (mmHg)\n\n","male_bp_abs.pdf",0.1,-0.93,10.5,"#B63679FF",-25,25)
sysdiff_male <-output%>%mutate(Male=sex)%>%select(-sex)
bwMsysdiff  <- bw
plotsysdiffM<-color_plot
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```

```{r ,results= TRUE, fig.height=22}



figuresysdiff <- ggarrange(plotsysdiffT, plotsysdiffF, plotsysdiffM,
                    labels = c("Total population", "Females", "Males"),
                    ncol = 1,
                    font.label = list(size=28, face="bold", family="serif"),
                    label.x = c(-0.12,-0.03,-0.01), label.y=c(1.08,1.08,1.08,-10,10)) +
                    
  theme(plot.margin = margin(2,1,2,1, "cm")) 

figuresysdiff
```

\newpage

**Outcome : Diastolic**

**Females and males**  

```{r tfull}

carrs_rdd(data$diasdiff,data$runningvar,data,"Change (mmHg)\n\n","overall_diastolic_bp_abs.pdf",0.16,-0.83,16.8,"#25AC82FF",-25,25)
diasdiff_total <-output%>%mutate(Total=sex)%>%select(-sex)
plotdiasdiffT<-color_plot
bwTdiasdiff  <- bw
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```

\newpage

**Outcome : Diastolic**

**Female**  

```{R tfemale}
data1<-data[female==1,]

carrs_rdd(data1$diasdiff,data1$runningvar_female,data1,"Change (mmHg)\n\n","female_diastolic_bp_abs.pdf",0.16,-0.91,16.8,"#46337EFF",-25,25)
diasdiff_female <-output%>%mutate(Female=sex)%>%select(-sex)
plotdiasdiffF<-color_plot
bwFdiasdiff  <- bw
rm(bw)

```

```{r ,results= TRUE}
summary(out)
```

\newpage

**Outcome : Diastolic (Male)**

```{r tmale}
data2<-data[female==0,]

carrs_rdd(data2$diasdiff,data2$runningvar_male,data2,"Change (mmHg)\n\n","male_diastolic_bp_abs.pdf",0.16,-0.95,16.8,"#B63679FF",-25,25)
diasdiff_male <-output%>%mutate(Male=sex)%>%select(-sex)
plotdiasdiffM<-color_plot
bwMdiasdiff  <- bw 
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```

```{r ,results= TRUE, fig.height=22}
figurediasdiff <- ggarrange(plotdiasdiffT, plotdiasdiffF, plotdiasdiffM,
                    labels = c("Total population", "Females", "Males",0,150),
                    ncol = 1,
                    font.label = list(size=28, face="bold", family="serif"),
                    label.x = c(-0.12,-0.03,-0.01), label.y = c(1.1,1.1,1.1)) +
                    theme(plot.margin = margin(2,1,2,1, "cm")) 
```


```{r ,results= TRUE, fig.height=25,fig.width=20}
figurecomb <- ggarrange(figuresysdiff, figurediasdiff,
                    labels = c("Systolic", "Diastolic"),
                    ncol = 2,
                    font.label = list(size=32, face="bold", family="serif"),
                    label.x = c(0.4,0.37), label.y = c(1.03,1.03)) +
                    theme(plot.margin = margin(2,1,2,1, "cm")) 
figurecomb
filename <- paste0("combined_bpdiff.pdf")
    ggsave(filename, dpi = 300)
    
filename <- paste0("combined_bpdiff.jpeg")
    ggsave(filename, dpi = 300)

```


```{r}

total_table<-merge(sysdiff_total,diasdiff_total, by="rowname")
total_table<-total_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(total_table)[2] <- "Systolic (change)"   
colnames(total_table)[3] <- "Diastolic (change)"   

female_table<-merge(sysdiff_female,diasdiff_female, by="rowname")
female_table<-female_table%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(female_table)[2] <- "Systolic (change)"   
colnames(female_table)[3] <- "Diastolic (change)"  

male_table<-merge(sysdiff_male,diasdiff_male, by="rowname")
male_table<-male_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(male_table)[2] <- "Systolic (change)"   
colnames(male_table)[3] <- "Diastolic (change)"  
```
##table making function

```{r}
#combining Systolic and Diastolic
  filename <- paste0("mainresults_combined_bpdiff.tex")
  caption=paste0("Effect of home-based screening on changes in blood pressure")

  comb_table<-rbind(total_table,female_table,male_table)

  combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccc", format = "latex",col.names = c(" ","Systolic","Diastolic"),digits = 2) %>%
      kable_styling(latex_options = "HOLD_position", position = "center") %>%
      pack_rows("Panel A: Total", 1,5) %>%
      pack_rows("Estimation results", 1,3) %>%
      pack_rows("Bandwidth details", 4, 5) %>%
      pack_rows(" ", 6,10) %>%
      pack_rows("Panel B: Female", 6,10) %>%
      pack_rows("Estimation results", 6, 8) %>%
      pack_rows("Bandwidth details", 9, 10) %>%
      pack_rows(" ", 11,15) %>%
      pack_rows("Panel C: Male", 11,15) %>%
      pack_rows("Estimation results", 11, 13) %>%
      pack_rows("Bandwidth details", 14, 15) %>%
      footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average Diastolic effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
      save_kable(filename, format = "latex")

```
  

